// dtrt.C -- Extract the drift times, and integrate to obtain the rt relations
// for dcs around the target
// Created: Jul 7, 2024 by SUN Yazhou

void a(vector<TH1 *> &v, TH1 *h){ v.push_back((TH1 *)h); }
void dtrt(){
  TFile *f = new TFile(getenv("TMP1"));
  TTree *etf = (TTree *)f->Get("etf");
  vector<TH1 *> v, vrt;
  TH1 *h0, *h1;
  char zcut[64] = "";
  static const int Z = 29;
  if(Z > 0) sprintf(zcut, "az0>%.2f && az0<%.2f", Z-0.5, Z+0.5);

  const char ud[2] = {'U', 'D'}, xy[2] = {'X', 'Y'}, *trkr[] = {"pre", "post"};
  char hname[64]{}, layer[64]{};
  const double RG[2][2] = {{-50., 200.}, {-100., 400.}}, MAXR[2] = {2.5, 5.};
  for(int i = 0; i < 2; i++){ // loop over U, D
    for(int j = 0; j < 2; j++){ // loop over 0, 1
      sprintf(layer, "dcTa%c%d", ud[i], j);
      sprintf(hname, "h%s", layer);
      h0 = new TH1F(hname, Form("Drift Time for %s", layer), 2000, RG[i][0], RG[i][1]);
      h0->GetXaxis()->SetTitle("drift time [ns]");
      h0->GetYaxis()->SetTitle("counts");
      h1 = new TH1F(Form("hrt_%s", layer), Form("R(t) for %s", layer), 2000, RG[i][0], RG[i][1]);
      h1->SetLineColor(2);
      h1->GetXaxis()->SetTitle("drift time [ns]");
      h1->GetYaxis()->SetTitle("drift distance [mm]");
      a(v, h0); a(vrt, h1);
      for(int k = 0; k < 2; k++) for(int l = 0; l < 2; l++)
        etf->Project(Form("+%s", hname), Form("%s%c%d.GetT()", layer, xy[k], l),
          // fired(): 2: normal; 3: global tracking
          Form("%s&&%s%c%d.fired()>=2&&fabs(%sTaTracker%c.ang())<1.", zcut, layer, xy[k], l, trkr[i], xy[k]));
      for(int i = 1; i < h0->GetNbinsX(); i++) h1->SetBinContent(i, h0->Integral(1, i));
      h1->Scale(MAXR[i]/h0->Integral());
      cout << "Processing " << layer << ": " << zcut << "...\r" << flush;
    } // end for over 1, 2
  } // end for over U, D
  cout << endl;

  TCanvas *c = new TCanvas("c", "Integrate t for R(t) for dcTa-s", 1680, 1050);
  c->Divide(2, 2);
  int id = 0;
  c->Print("dtrt.pdf[");
  for(int i = 0; i < 2; i++){ // U,D
    for(int j = 0; j < 2; j++){ // 0, 1
      c->cd(2*j+1); v[id]->Draw();
      c->cd(2*j+2); vrt[id++]->Draw();
    } // end for over j
    c->Print("dtrt.pdf");
  } // end for over i
  c->Print("dtrt.pdf]");

  TFile *frt = new TFile("rt.root", "UPDATE");
  for(auto h : vrt) h->Write("", TObject::kOverwrite);
} // end dtrt
